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ABSTRACT 

An unsteady compressible viscous waJte flow past a circular cylinder has been successfully 
simulated using spectral methods. A new approach in using the Chebyshev collocation method for 
periodic problems is introduced. We have further proved that the eigenvalues associated with the 
differentiation matrix are purely imaginary, reflecting the periodicity of the problem. It had been 
shown that the solution of a model problem has exponential growth in time if “improper” boundary 
conditions are used. A characteristic boundary condition, which is based on the characteristics of 
the Euler equations of gas dynamics, has been derived for the spectral code. The primary vortex 
shedding frequency computed agrees well with the results in the literature for Mach = 0.4 , Re = 80. 
No secondary frequency is observed in the power spectrum analysis of the pressure data. 


‘Research was supported by the National Aeronautics and Space Administration under NASA Contract Nos. 
NASl-18107 and NASl-18605 while the authors were in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1 Introduction 

The behavior of the flow past a circular cylinder has been used as the model for fundamental studies 
of external (open) flows (figure 1.1 shows the schematic of the flow past a circular cyUnder). In 
spite of the simplicity of the geometry, it embodies many interesting features of fluid dynamics, 
ranging from steady Stokes flow, viscous wake flow to turbulent flow [1]. Morkovin [2] described 
the behavior of flow past a circular cylinder as the “kaleidoscope of challenging fluid phenomena”. 

Over the years, the unsteady viscous compressible wake flow past a circular cylinder has been 
extensively investigated both theoretically and experimentally by many researchers. Theodore 
von Kaxmtin [3] showed that the formation of the inviscid Karman vortex trail at low Reynolds 
number (see section 2 for definition) behind a bluff body moving downstream with the flow is 
linearly stable. The primary vortex shedding frequency is well established both experimentally 
and numerically over a wide range of Reynolds numbers [4]. If the frequency is nondimensionalized 
by the free-stream velocity and the diameter of the cylinder, it is called the Strouhal number St 
in the literature. For example, the Strouhal number of primary vortex shedding frequency Sti at 
Reynolds number Re=80 is about 0.156. 

A recent paper by Sreenivasan [5] reported experiments indicating that the primary shedding 
frequency is modulated by a lower secondary frequency, which is not a sub-harmonic to the shedding 
frequency (see figure 1.2). He suggested that this lower frequency played a major role in the 
period doubling route to chaos. Moreover, he also observed a discontinuity of the vortex shedding 
frequency versus Reynolds number (windows of chaos), which was observed earlier by Tritton [6, 7], 
Berger and Wille [8] and Priehe [9]. Sirovich [10] suggested some possible mechanisms to justify the 
existence of the secondary frequency. He pointed out that the perturbation and interaction among 
the vortex pairs can create a secondary frequency. However, Sreenivasan *s finding was disputed by 
Van Atta and Gharib [11], who attributed it to aeroelastic coupling of the flow with the cylinder’s 

vibrating modes. 

In light of all these contradictory arguments of the secondary frequency among experimental- 
ists, numerical simulations seem to offer an alternative means of gaining some insight into this 
controversy. The numerical simulation of unsteady compressible wake flow past a circular cylin- 
der by Townsend, Rudy and Sirovich [12] using a second-order finite-difference scheme seems to 
confirm the existence of this controversial secondary frequency. The time trace of the pressure 
and the power spectrum analysis (figure 1.3) clearly shows that the primary shedding frequency 
was modulated by a lower secondary frequency. The numerical value of the secondary frequency 
computed, however, is domain size dependent. Karniadakis and Triantafyllou [13] simulated the 
incompressible flow past a cylinder using the spectral-element technique. No secondary frequency 
was observed in their simulation. A secondary mode was excited by introducing the forcing term 
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into the momentum equation. 

In this paper we report the results of a collocation (pseudospectral) simulation of the compress- 
ible viscous Navier-Stokes equations. Our results demonstrate clearly that the secondary frequency 
does not exist. However, we had diihculties in imposing outflow boundary conditions in a stable 
way. In particular the boundary conditions used in [12] yield instabilities in the spectral code. 
The only set of boundary conditions that stabilized the code will be described in detail and model 
problems will be analyzed. 

When these boundary conditions were implemented in the finite-difference code as well the 
secondary frequency disappeared. The secondary frequency found in the finite-difference code 
had been a purely numeric^ phenomenon. Moreover, Sreenivasan reported that changing the 
aspect ratio of his wind tunnel, in effect pushing the outer boundaries further away, eliminated the 
secondary frequency in the experiment. 

This research illustrates an important feature of spectral methods. In contrast to the low 
order (and more robust) finite-difference schemes, spectral methods are less tolerant to incorrect 
implementation of boundary conditions. Thus the same boundary conditions that produced stable 
but erroneous results in the finite-difference code were not tolerated in the spectral code. Results 
with spectral codes are harder to achieve but they are reliable! 

This paper is structured as follows: Section 2 describes the governing equations. In section 3 
we describe the spectral code based on Chebyshev and Fourier collocation methods. In this section 
we describe also the filtering used and several options of grid generation. ATs^in this section we 
present a proof that the eigenvalues of the Chebyshev method applied to a periodic problem are 
purely imaginary. Section 4 discusses the boundary conditions, this section is really the heart of 
this paper. In section 6 we present the numerical results. The conclusion is presented in section 6. 


2 Navier-Stokes Equations 


The governing equations are the two-dimensional compressible viscous Navier-Stokes equations in 
strong conservation form [14]. They describe the laws of conservation of mass, momentum and 
energy in the absence of external forces. This set of equations in the Cartesian coordinates {x , y) 
is transformed into the general curvilinear coordinates , r}) with the transformation Jacobian 
^ — isVy ~ The nondimensional form of these equations is 


dF dG _ 1 (dPu doA 
dt ^ ^ dr] ~ Re \d^ ^ dri ) 


( 2 . 1 ) 
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coupled with the nondimensional equation of state and the internal energy, respectively, 

p = ( 7 -l)pT, E = p]^T + ^{u^+v^)^ . 

The elenaents of the stress tensor are 
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and 


dT _dT^ ^£2 ^ = + 

aT ~ "d^dx ^ drjdx’ dy dC dy ^ drj dy ’ 

The viscosity is related to the temperature by the Sutherland viscosity law 


M = 


( 2 . 2 ) 

Cj + T 

where and c, are constants. The ratio of specific heats 7 is 1.4 and the Prandtl number Pr is 
0.72. Cy is the specific heat at constant volume and k is the coefficient of thermal conductivity. 

The free-stream Reynolds number Re = p^oUcoDlfioo of the flow is based on the diameter of 
the cyUnder D, free-stream velocity Uoo, free-atream density poo and dynamic viscosity poo- The 
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parameter Re in the Navier-Stokes equations (2.1) is the reference Reynolds number computed 
by the reference quantities Cf„/ = C^oo , Pref = P«, , ?;«/ = £^i/cv • We shall refer to the 

free— stream Reynolds number as Re unless clarification is deemed necessary. 

The Navier-Stokes equations (2.1) and their necessary supplementary relations such as the 
Sutherland viscosity law (2.2) are characterized by four parameters [14]. These parameters are the 
free-stream mach number Moo, free-stream Reynolds number Re, diameter of the cylinder D, and 
dimensional temperature Tq. 

3 The Numerical Method 

The 2-D Navier-Stokes equations are solved numerically in polar coordinates with angular (0) 
and radial (r) coordinates. Instead of using mapping techniques, such as exponential or algebraic 
transformations [16, pages 71-76], to map the infinite physical domain into a finite size, the physical 
domain is simply truncated at some fixed distance away from the cylinder center (see figure 3.1). 
Typically the outer computational boundary is located about 23 cylinder diameters away from the 
center of the cylinder. 

An 0-type grid configuration with branch cut in the center of the wake behind the cylinder is 
mapped into the computational domain with coordinates f = ^(0) and r} = vi’’’), j'6-i 

X = r{ri) cos d(f) , y = r(»/) sin 0{$) . 

The cylinder wall boundary, the outer computational boundary and the periodic-cut boundary 
are mapped into the four boundaries of the computational domain. This choice of grid structure 
offers two major numerical advantages. First, the implementation of the boundary conditions is 
simplified. Secondly, this choice of coordinate system is particularly suited for the employment of 
the Fourier collocation method in the angular periodic direction and the Chebyshev collocation 
method in the radial non-periodic direction. They are employed in the form of matrix vector 
multiplication method instead of the commonly used Fast Fourier Transform method. The reader 
is referred to [15, 16, 17] et al. for a detailed discussion of these methods. 

In some of the numerical experiments, we employed a commonly used exponential filter 

r 1 ^ 0 < |Jfc| < ko 

“ I ko < |ife| < 

where kc is the cut-off frequency, to improve the stability of the spectral scheme. The parameter 
a is chosen in such a way that it drives to machine zero e, i.e., o = - Ine. The parameter 7 
(7 “ 4 or 6 in this research) determines the order of the filter. The filter is incorporated into the 
solution and differentiation matrices for both the Fourier and Chebyshev collocation methods. 

In this application, the important physical phenomena are confined to the viscous wake behind 
the cylinder and the boundary layer on the cylinder surface. In order to improve the resolution 
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within these regions without resorting to allocating additional grid points, the grid transforma- 
tion technique is used for grid point redistribution. Since we found that Chebyshev collocation 
points provide sufficient resolution within the boundary layer, no additional mapping technique is 
employed in the radial direction. 

The angular grid transformation [18] used is expressed in the complex form 


je _ 




(3.2) 


where /3 is a parameter controlling the amount of clustering of grid points at location 7 G [0,27r], 
and |/3| < 1. 

The advantages of this grid transformation are summarized as follows : 


1. 9{^) and all of its derivatives are analytic and periodic. Hence, the periodic boundary condi- 
tions are satisfied and spectral accuracy is retained for spectral methods. 

2. 0(^) preserves the periodicity of the Navier-Stokes equations in the angular direction, 

3. Analytical transformation derivatives are readily available without resort to numerical ap- 
proximation, which is usually done in grid generation techniques. 


In this numerical simulation of flow past a circular cylinder, /? = —0.4 with 7 = 0 provides 
a satisfactory resolution within various flow regions in the angular direction. At least half of the 
grid points used are clustered within the wake region. The other half are distributed around the 
remaining regions with gradually reduced resolution toward the front of the cylinder (see figure 
3.1). 

A seemingly attractive alternative to using the Fourier collocation method with mapping is 
to use the Chebyshev collocation method. As noted before the Chebyshev collocation points are 
unevenly distributed and cluster at both end points with spacing proportional to 1 / , In our case 
the two end points coalesce to 0 = 0 and this is exactly the region where we need more resolution. 

Here we would like to further investigate some aspects of the application of Chebyshev methods 
for periodic problems. There is more than one way of implementing the above idea. We will discuss 
one method used for discretization of the first derivative operator and show that it leads to a 
differentiation matrix with purely imaginary eigenvalues, reflecting the periodicity of the problem. 

Consider the model problem 


Ut 

= Ux 


= 


= /(*)• 


(3.3a) 

(3.3b) 

(3.3c) 


The simplest way to apply the Chebyshev collocation (pseudospectral) method to the model 
problem is to satisfy the equation (3.3a) at the points {x^ = cos {'Kj/N) , j = 1, . . . , i\T - 1} 
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and to impose the periodicity condition (3.3b). Since the numerical approximation Uni^x^i) is a 
polynomial of degree N , we need another condition, usually taken to be 


To summarize Un satisfies 


dUtf 

dx 


i-ht) 


dx 


(1.0 


dUff 

dt 


9uif 

dx 


+ (Ax + B)Tj, 


Urr{-l,t) = itAr(l,0 


(3.4) 


(3.5a) 

(3.5b) 


For even N, is an even function and is odd, therefore, 5 = 0 in (3.6a). 

We are ready now to state 

Theorem 1 The eigenvalues of the periodic Chebyshev approximation to the first derivative matrix 
are purely imaginary. 


Proof : Let y? (x, a) be an eigenfunction of the first derivative periodic operator and s an eigenvalue 
then 


dip 

dx 

¥’(-1,0 


sip-xTlf 
¥>(l,s) . 


(3.6a) 

(3.6b) 


Since y>(x,s) is a polynomial of degree i\T in x, it can be written explicitly as 


¥’(®.0 = Z ) 


Jfc=0 


,k+l • 


(3.7) 


We denote by 


Thus 


^(i/) = vip{l, -) 


/c=0 


*=0 


k=0 


and finally upon setting ^ we get 


V'(J/) = h{fi) + vg{/i) 


(3.8) 
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where 


g(^) = 

Jb=0 

k(ii) == f;(.r„)<“)(i)/ 

ife=0 . 

It had been shown that the roots of V* ( »^ ) have all negative real part, 
of the stability of the non-periodic Chebyshev method. By the theory of 
roots of g{fi) and h ( /x ) are real and negative. 

Substituting (3.7) into (3.6b) we get that the eigenvalues 3 of the 
Chebyshev polynomials satisfy 

k=0 ^ 

Now since sc is an even polynomial only the even powers of s remain in 
(3,9) can be rewritten as = 0. The roots fi are real and negative, 

s are purely imaginary. This concludes the proof. 

A very surprising result follows from a closer look at h(/i) = 0. This equation is the same 
as the characteristic equation for the Chebyshev second derivative matrix with Dirichlet boundary 
conditions. Thus the eigenvalues of the first derivative periodic Chebyshev are the square roots of 
those of the second derivative matrix with Dirichlet boundary conditions. 

We plan, in the future, to further investigate the use of the periodic Chebyshev method for the 
azimuthal direction. In this research we used standard Chebyshev in the radial direction and the 
mapped Fourier method in the azimuthal direction. 

4 Boundary Conditions 

In this section, we examine the procedures of imposing various boundary conditions on the outer 
computational boundary and the cylinder surface. Due to the geometry of this application, the 
Fourier collocation method enforces the periodicity of all variables automatically in the angular 
direction. The only challenge is to find a stable procedure in providing boundary conditions at the 
cylinder surface and the outer computational boundary in the radial direction. 

Since, in the present research, the viscous regions are concentrated within the thin boundary 
layer on the cylinder surface and the narrow near walce behind the cylinder, it is quite reasonable 
to assume that the flow is basically inviscid at the outer boundary if it is located sufficiently far 
away from the cylinder. At the cylinder surface, however, viscous flow boundary conditions are 
used. Also for the present research, the flow is subsonic everywhere. Thus the boundary conditions 
can be taken from the analysis of hyperbolic equations. 


this being a consequence 
Hurwitz polynomials the 

first derivative periodic 
(3.9) 

(3.9). Denoting fi = 
therefore the eigenvalues 
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It is well known that great care should be exercised in the application of boundary conditions 
when spectral methods are used for numerically solving hyperbolic systems of equations. We wish 
to review here the theory of Chebyshev collocation methods applied to hyperbolic systems. Rather 
than quote the theorems we wish to illustrate the theory by examples. 

Example 1 : Consider the system of equations 


Ut = Ua. + Ux 

Vt = - Wa, (4.1) 

with the boundary conditions 


= 0 


(4.2) 


The straightforward application of the Chebyshev method can be summarized in three steps ; 

A ) Given u{xj,t),v{xj,t) Xj = cob^-kJ/N) j = 0, . . . ,iV form and |^(aj,-,t) . 

B ) Use a time marching technique to get u(iy, t + At) and v{xj,t + At) again for j = 0, . . . , JV . 
C ) Change u{xo,t + At) and u{xn,t + At) to satisfy the boundary conditions (4.2). 

Thus the values of v at the boundaries are taken from step B. 

We wish to demonstrate that the procedure outlined above is unstable. Consider the eigenvalue 
associated with (4.1) : 


sil = Ua: + Vx + {Ax + B)T^ 
sv = -fix 
fi(±l,s) = 0 


(4.3a) 

(4.3b) 

(4.3c) 


where u and fi are i\Tth degree polynomials in x. Equation (4.3a, b,c) reflects the fcict that the 
boundary conditions are imposed on the variable u only. 

Close inspection of (4.3b) reveals that fi = 0. The reason is that the left hand side of (4.3b) is 
a polynomial of degree N whereas the right hand side is a polynomial of degree N - 1. Thus we 
are left with the equation 


sil = Ux 4- fix + {Ax + B)T'f, 
fi(il,3) = 0 

A polynomial solution of (4.4a) can be expressed in the form 


OO 

ii{x,s) = 

Jb=0 


(xn)W 


OO 


+ SE 

fc =0 


^fc+1 


(4.4a) 

(4.4b) 
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Equation (4.4b) reads now 






teo 

“ (xn)W(-i) 


cfc +1 


k-0 
oo 


+ ^ 2 ^ jfc+i 


= 0 

(4.5a) 

= 0 

(4.5b) 


Upon adding and subtracting (4.5a) and (4.5b) and making use of the fact that x T'^f is even 
and T'ff is odd, one gets 


.4 E (» n)<“>(i) + B E = 0 

fc =0 k =0 


k=0 


fc =0 


where u = 1/s, 

We denote by 

Jb=0 

= EW)‘“+'’(i)‘''“> 

fc =0 

So that (4.6) can be written as 


k=0 


*=0 


A h^{u^) + = 0 
Avg,{v^) + B h,{u^) = 0 


(4.6) 


(4.7) 


(4.8) 


Thus the eigenvalues are determined by the characteristic equation 

Equation (4.9) is expressed in term of — /x. Thus the only possibility for eigenvalues s without 
positive real part is that the polynomial 

hi(/i)hj(/i) - MPi(/i)s,(/i) = 0 (4.10) 


has real negative roots. 

Equation (4.10) is a special case of equation (4.10) in [19] with q = /3 = 7 = ^. The solution in 
these cases must have positive real part since it approximates a solution of the form e . Thus v 
must have positive real part. This conclude the proof. 

Example 2 : Consider again equation (4.1)-(4.2). Define 

w{x,t) = 2u{x,t) + v{x,t) (4.11) 
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then it is clear that 


dw _ dw 
It ~ 

dv dv 

a = (^-12) 

with the boundary conditions 

u;(l,t) = v(l,t) 

= v(-l,<) (4.13) 

The problem (4.12)-(4.13) is identical with (4.1)-(4.2). However the application of the Cheby- 
shev collocation method is different than the one outlined before. Basically the first two steps of 
the algorithm described above remain unchanged. 

In the third step we change not only u{x,t) but also w(x,t) by using equations (4.11)- -(4.12) 
in the following way : 

At ® = — 1 one has to specify i;(-l,t) and to use «;(— l,t) as computed from steps A) and B). 
Thus, denoting by u;c(-l,t -j- At) the value of 2u(-l,t -h At) -b u(-l,t -j- At) we have the system 

2«(-l,t + At) -1- v(-l,t -}- At) = u;(-l,f -1- At) = 4- At) 

u(-l,t-|- At) = 0 (4.14) 

Equation (4.14) is a system of equations used to update ti(-l,t-|- At) ^ v(-l,t -|- At). A similar 
procedure is carried out at the point x = 1. 

The procedure outlined in (4.14) cast the numerical procedure in the form discussed in [20]. 
Of course the best procedure is to specify u;(l,t) and u(-l,t) and to use u;c(-l,t) and Uc(l,t). 

This procedure wiU be used in the current paper and referred to as the characteristic boundary 
conditions. 

4.1 Characteristic Boundary Conditions 

In this section, we will give the general procedure of how to impose the characteristic boundary 
conditions at the outer computational boundary followed by the specific implementation for the 
numerical experiment. 


General Procedure 


Consider the system of one-dimensional nonlinear hyperbolic partial differential equations with N 
dependent variables in the conservation form 


— + = 0 


(4.1.1) 


r i 
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and the Jacobian of F is defined as A{W) = . 

The procedure consists of three main steps : linearization, characteristic decomposition and 
characteristic treatment. 

Linearization : 

Equation (4.1.1) is linearized around some uniform state i-e., 


dW 

dt 


+ A{vTo) 


dx 


= 0 


(4.1.2) 


where the linearized Jacobian A(ulo) is a. constant N x N matrix. 

Characteristic decomposition : 

The system (4.1.2), which is a hyperbolic system with constant coefficients, is diagonalized into 
the characteristic form, i.e.. 


dRi ^ , dRi 
dt * dx 


= 0 , 


i= 1,2,. ..,N 


(4.1.3) 


where are eigenvalues and Ri = Ri(^ , t?o) are characteristic variables (eigenfunctions). 

The sign of each eigenvalue A» is then examined at each outer computational boundary point. 
Suppose that 

f A» > 0 i = 1, . . . , A (4 14) 

\ Ai < 0 i = k + ’ 

then k quantities must be specified and N — k quantities should be obtained numerically at that 
boundary point. 

Characteristic treatment : 

Let denote the numerical value at the outer computational boundary. Spectral methods 
update every grid point in the domain, including all of the boundary points, at each time step. 
No additional extrapolation procedure is needed for the boundary values. Furthermore, we denote 
u^a/ to be the known value of the far— field conditions, for example the free-stream conditions at 
infinity. 

According to the linear theory, the characteristic quantities .Ri(W^,n?o),» = l>’--)^ should be 
specified and the functions R{ ( ^ )i * — k 1, . . . ,N should be obtained from the numerical 
scheme. Thus 


Ri{\^,v}o) = Ri{^nu,^o), i = k + l N, (4.1.5) 

This process results in a linear system of N equations with N unknowns. Once solved, this 
system of equations yields W at each boundary point. We will henceforth denote this procedure 
as characteristic boundary conditions. 
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Characteristic Boundary Conditions for the 2-D Euler Equations of Gas Dynamics 

Following the general procedure outlined in the previous section, we will derive the characteristic 
boundary conditions for the 2~D Euler equations of gas dynamics. 

Let us denote 


^oo = («oo ,Voo) = (1 ,0) 

= (mu,mv) = {f>u,pv) 

c = y/lin - 1)^00 

^ = (i.i) 

= {vx ,vv) ! 


to be the free-stream velocity, 
to be the mass flux, 

to be the local free-stream sound speed, 


to be the unit radial normal vector, 
and (/>, TTiu, TOv, E) to be conserved mass, momentums and internal energy. The subscript oo 
denotes free-stream values. 

Linearization : 

Consider the 2-D Navier-Stokes equations (2.1), 


dW 

dt 


+ 


9F 

dr] Re \ dt] J 


where — (p, ttIu, ttIv, E)*. 

Far away from the boundary layer on the cylinder surface, the flow is dominated by the con- 
vective term and viscosity exerts only a minor effect except inside the narrow wake. Hence, by 
assuming u is small or Re is large, (2.1) becomes the 2— D nonlinear Euler equations in strong 
conservation form. The nonlinear Euler equations can then be linearized with the free-streaim 
conditions u/o = (poo > ^oo , -Eoo), i-e., 


dW , dW „ dW ^ 


(4.1.6) 


where Aqq and are constant Jacobian matrices in the radial and angular directions respectively. 
Since we are only interested in deriving the boundary conditions in the radial direction, we can 
drop the angular dependent term Bqq'^W from consideration. 

Characteristic decomposition : 

It is a matter of algebra to find the eigenvalues and their corresponding characteristic yari- 
ables Ri of .4oo* They are 


Ai = Aq = ifoo ' ^ ) A3 = if 00 * N — c y A4 = if 00 • + c , 


(4.1.7) 
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and 


Jii = j - Qpt/^oo-^oo - ^-^oo + 

iZa = (’’^U -P'“oo)|^ - (’’iv -pVoo)^ (4.1.8) 

Jt4 = + 

Characteristic boundary conditions : 

The values of all conservative variables are assumed to be known at the outer computational 
boundary. 

At the subsonic inflow region, i.e., ■^oo • y < 0 and ^oo • 1 ^ < c, one gets 

Ai<0, i = l,2,3 ; A4>0 

from (4.1.7). Accordingly, 

iij (p , ii? , E) — F% = R% {poo ) oo ) ^oo) j i = 1, 2, 3 

R4{P,H,E) = Fi = Ra {pnu , j^n« , (4-1-9) 

should be used at each inflow boundary point. 

At the subsonic outflow region, i.e., ~lfoo • ^ >0 and ^oo ■ ^ ^ gets 

A*>0, i= 1,2,4 ; A 3 < 0 

from (4.1.7) instead. Accordingly, 

R{{p^lliiyE) = Fi R{ {Pnu J nu i ^nu) i i= 1,2,4 

R 3 {p,Ji,E) = Fi R 3 (Poo , ^00 , E^) (4.1.10) 

should be used at each outflow boundary point. 

Remark : Since the grid structure can be constructed in such a way that if (0 , 1) at any outer 
boundary point, the case of if^o -if = 0 mil never arise in this application. 

More specifically, equations (4.1.8), (4.1.9) and (4.1.10) together yield 
Fx = p^-{^plfco-lfoo-^-lfoo + E^ 
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^2 = 

(ro„ - /) o„) - (m, - )j ^ 

(4.1.11) 

Fz = 

- (S - , I7„) . V + C' - 1) (i„ ZL -tt-t.+E) 


Fa = 

{J^-pVoo) ■ C^oo - + E^ . 



It is easy to check that if 

= ,f>, = l{F,-F3) , h = F, , 

then the solution of (4.1.11) is 

P = 

m„ = (^3% + ^3^) +/>««, (4.1.12) 

m„ = {<f>2%- pvoo 
E = <j>i- -plfoo-Tfoo + m„«oo + m„Voo . 

The primitive variables {p,u,v,T) can then be extracted from (4.1.12) at each outer computational 
boundary point. 

4.2 Cylinder Surface Boundary Conditions 

At the cylinder surface, the no— slip and constant temperature boundary conditions are enforced, 
i.e., 

u = u = 0, T = Too, P = Pc^, (4.2.1) 

where p^^ is the density computed numerically by the spectral code. 

5 Numerical Results and Discussion 

The spectral code is fully vectorized to taie advantage of the vector pipeline hardware architecture 
of the Cyber 205 supercomputer at the NASA Langley Research Center. Baaed on the free-stream 
flow conditions {Uoo = 186.0 m/sec. Too = 538.3 IT, Poo = 30.61 Pa) and the cylinder diameter 
D = 61.0 mm, all the computations are made at free-stream Mach number, Mach = 0.4, and 
free-stream Reynolds number, Re = 80. 

To advance the flow field in time, the second-order Runge-Kunta method is employed. The 
initial flow is started impulsively (uniform free-stream conditions throughout the whole domain 
except the no— slip and constant temperature boundary conditions at the cylinder surface). After 
the transient efifect of the initial conditions has diminished from the flow— field, the numerical 
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scheme is integrated over 20,000 or more time steps providing enough cycles of vortex shedding for 
analysis. The time trace of the pressure signal at several locations is recorded and stored for later 
power density spectrum analysis using the Discrete Fast Fourier Transform (DFFT) technique. 
The chosen point for which the signal was analyzed is located at approximately ID above the wake 
centerUne and lOD behind the cyHnder within the wake (X,0- We also recorded the pressure 
signal at location lOD away from the cylinder center and 180 degrees around the cylinder in the 
inflow re^on (Xi„). In order to be able to compare the numerical results with the experimental 
results, the signal point X*k is located approximately at the position of the hot wire probe in some 
experiments of Sreenivasan. 

Each set of time series contains the nondimensionalized pressure P{t)l Foo s.t each time step. 
Since we are only interested in the perturbation signal superimposed on the uniform flow, the 
average of P{t)fPoo is first subtracted from each data set. The resulting time series P{t) is analyzed 
in the frequency domain in order to identify individual dominant frequencies. This is usually 
accomplished by the power density spectrum analysis. A discrete time series P{tn) of length T is 
transformed into the frequency space by the Discrete Fast Fourier Transform technique. 

If /, A: and At denote the frequency, the wave number and the time increment respectively, 

then 

PM = 

k=-N 

where T = NAt, tn = nAi and k = fT. 

For each wave number k, the complex coefficient a* is multiplied by its conjugate 3*, i.e. 

= 0 * 3 *. {6*, k = 0,. .. ,N} is called the power spectrum, which measures the energy of each 
wave number (frequency). Normally, the power of each wave number is plotted using a logarithmic 
scale of base 10. Those dominant frequency components will show up as sharp peaks in the plotted 
power spectrum. Since the DFFT of the pressure signal requires a uniform time increment At, it 
is fixed at 0.00434201 after the flow is settled in a regular pattern of alternate vortex shedding. 

Spectral Scheme with characteristic boundary conditions - SPl 

The resolution of the spectral code, which used 64 X 48 grid points in the angular and radial 
directions respectively, is comparable to, if not better than, the resolution of the full finite— difference 
code [12], which used 120x 151 gird points. The spectral scheme becomes unstable when the number 
of grid points used is reduced in either directions. It appears that 64 x 48 is a lower limit of the 
number of grid points required to achieve the resolution. This conforms to the typical exponential 
convergence behavior of spectral methods. This should come as no surprise because the spectral 
code is developed within the theoretical scope of spectral methods. Furthermore, since a lesser 
number of grid points is employed in the spectral scheme, the larger spacing also allows a larger 
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time step to be taken in the numerical simulation. Typically, the time step of the spectral code is 
about 10 times larger than the full finite-difference method. 

Figures 6.1a-b show contour plots of density p/poo and temperature TfT^, respectively. The 
vortices pass through the outer computational boundary without any apparent reflection back 
into the interior domain. Furthermore, no observable spurious perturbation travels back into the 
interior domain either. The Karman vortex street is clearly demonstrated in the contour plots. 

The pressure time series and its power spectrum at X„k and Xi„ are shown in figures 5.2a-b 
respectively. As indicated in all the power density spectra, no secondary frequency can be observed 
both within the wake and in front of the cylinder at the inflow region. The frequencies that appear 
in the power spectrum are the primary vortex shedding frequency Stj = 0.1574 and its harmonics 
only. 

Spectral Scheme with non— characteristic boundary conditions — SP2 

To show the sensitivity of the boundary conditions for spectral methods, we replace the character- 
istic boundary conditions by the non-characteristic one used in [12]. The initial flow is restarted 
from the flow— field data of the previous run SPl. The simulation becomes explosively unstable 
after only 71 time steps. 


Scheme 

Domain 
size (D) 

Number of 
time steps 

Total 

time 

Sti 

St2 

SPl 

46 D 

20480 

88.92 

0.1574 

none 

SP2 

46 D 

71 

unstable 



SPl 

41 D 

19456 

70.26 

0.1565 

none 


Table 1: Results of the spectral simulation of unsteady compressible flow past a circular cylinder 
with Nfl = 64, Nr = 48, Mach = 0.4, Re = 80. 


6 Conclusion 

An unsteady compressible viscous wake flow past a circular cylinder has been successfully simulated 
using the Fourier and Chebyshev collocation methods. A new approach in using the Chebyshev 
collocation method for periodic problems is introduced. We have further proved that the eigenvalues 
associated with the differentiation matrix are purely imaginary, reflecting the periodicity of the 
problem. It had been shown that the solution of a model problem has exponential growth in 
time if “improper” boundary conditions procedure is used. Hence, we derived the characteristic 
boundary conditions, which is based on the characteristics of the Euler equations of gas dynamics, 
for the spectral simulation. The primary vortex shedding frequency agrees well with the results in 
the literature for Mach number 0.4 and Reynolds number 80. No secondary frequency is observed 
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in the power spectrum an 2 ily 8 is of the pressure data that is generated by the spectral code with 
the characteristic boundary conditions. 
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Figure 1.1: Schematic of flow past a circular cylinder of diameter D. The location of the hot wire 
measurement in Sreenivasan’s experiments is marked as * . 



Figure 1.2; The power density spectrum of velocity u at Re=120, x/D - 8 as reported by Sreeni- 
vasan [ 5 ]. The inset shows a typical hot wire aerometer oscillograph. The power is plotted on a 
logarithmic scale of base 10. 
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Figure 1.3: The pressure time series and the power density spectrum of the full finite-difference 
scheme [ 12 ]. 



Figure 3.1: The grid structure of flow past a circular cylinder in Cartesian coordinates (i, y). 


20 



ORiGIN'AL PAGE IS 
OF POOR QUAUTY 



Figure 5.2a: Pressure time series and power density spectrum of the spectral scheme SPl at the 
point located 102? behind the cylinder and 12? above the wake centerline. 



sicr KHcn 


Figure 5.2b: Pressure time series and power density spectrum of the spectral scheme SPl at the 
point located 102? away from the cylinder center and 180 degrees around the cylinder in the inflow 
region. 
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